General relativistic simulations of magnetized binary neutron star mergers 
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Binary neutron stars (NSNS) are expected to be among the leading sources of gravitational 
waves observable by ground-based laser interferometers and may be the progenitors of short-hard 
gamma ray bursts. We present a series of general relativistic NSNS coalescence simulations both 
for unmagnetized and magnetized stars. We adopt quasiequilibrium initial data for circular, ir- 
rotational binaries constructed in the conformal thin-sandwich (CTS) framework. We adopt the 
BSSN formulation for evolving the metric and a high-resolution shock-capturing scheme to handle 
the magnetohydrodynamics. Our simulations of unmagnetized binaries agree with the results of 
Shibata, Taniguchi and Uryu In cases in which the mergers result in a prompt collapse to a 
black hole, we are able to use puncture gauge conditions to extend the evolution and determine the 
mass of the material that forms a disk. We find that the disk mass is less than 2% of the total 
mass in all cases studied. We then add a small poloidal magnetic field to the initial configurations 
and study the subsequent evolution. For cases in which the remnant is a hypermassive neutron 
star, we see measurable differences in both the amplitude and phase of the gravitational waveforms 
following the merger. For cases in which the remnant is a black hole surrounded by a disk, the disk 
mass and the gravitational waveforms are about the same as the unmagnetized cases. Magnetic 
fields substantially affect the long-term, secular evolution of a hypermassive neutron star (driving 
'delayed collapse') and an accretion disk around a nascent black hole. 

PACS numbers: 04.25.D-,04.25.dk,04.30.-w 



I. INTRODUCTION 

There is great interest in studying the effects of mag- 
netic fields in relativistic astrophysics. Magnetic fields 
are always present in astrophysical plasmas, which are 
usually highly conducting. Even if the initial magnetic 
field is small, it can be amplified via magnetic winding 
and other magnetic instabilities (see, e.g., @, Q for re- 
views). Neutron stars (NSs) have the strongest observed 
magnetic fields (up to ~ 10 15 G) among astrophysical ob- 
jects The strong magnetic fields result from gravi- 
tational collapse, which amplifies the magnetic fields in 
the core of the progenitor, and from various dynamo 
processes after the collapse (see, e.g., [f| for a review). 
Strong magnetic fields in a NS may trigger observable 
events such as pulsar glitches and the emission of large 
bursts of gamma rays and X-rays as in a soft gamma- 
ray repeater. Magnetic fields in a binary neutron star 
(NSNS) system may also influence the dynamics of the 
remnant after the NSs merge. 

Mergers of binary neutron stars (NSNSs) are expected 
to be among the leading sources of gravitational waves 
observable by ground-based laser interferometers. Ob- 
servations of short-hard gamma-ray bursts (GRBs) sug- 
gest that a substantial fraction of them may be formed 
from mergers of NSNSs or mergers of neutron star-black 
hole binaries (BHNSs). Many theoretical models of GRB 
engines consist of magnetized accretion disks around a 
spinning black hole [aS HI • General relativistic magne- 
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tohydrodynamics (GRMHD) is necessary to model such 
systems. 

The first two GRMHD codes capable of evolving the 
GRMHD equations in dynamical spacetimes were devel- 
oped by Duez et al. 0] (hereafter DLSS) and Shibata 
& Sekiguchi 

[3 (hereafter SS). These codes are based 
on the BSSN (Baumgarte-Shapiro-Shibata-Nakamura) 
scheme to integrate the Einstein field equations, a high- 
resolution shock-capturing (HRSC) scheme to integrate 
the MHD and induction equations, and a constrained 
transport scheme to enforce the "no-monopole" mag- 
netic constraint. Subsequently, Giacomazzo & Rez- 
zolla [ll| and Anderson et al. |12j developed similar 
codes. Our code (DLSS) and the code of SS have been 
used to study magnetic fields in hypermassive neutron 
stars [ill, 0, EBl > magnetorotational collapse of massive 
stellar cores to neutron stars llol . and (unmagnetized) 
coalescing BHNSs 0, H [IlTIoj]. We have also used 
our code to study the magnetorotational collapse of mas- 
sive stellar cores to black holes [21| as well as coalescing 
BHBHs [Hj], which are pure vacuum simulations. Shi- 
bata et al. have also performed simulations of (unmagne- 
tized) coalescing NSNSs [3, [H HI US HI] • 

Recently, Anderson et al. have used their code to study 
the coalescence of both unmagnetized and magnetized 
NSNSs [27], |28|. In the unmagnetized cases, they find 
an initial configuration in [27| that leads to prompt col- 
lapse to a black hole following the merger. The to- 
tal (baryon) rest mass of their initial configuration is 

M w 1.03M^ TOV) [H, where M^ TOV) is the maxi- 
mum rest mass of a nonrotating neutron star for the 
n=l polytropic equation of state (EOS) adopted in their 
simulation. This result seems to contradict the earlier 
finding of Shibata, Taniguchi & Uryu [l| that prompt 
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black hole formation is possible for this EOS only if 
M > 1.7M^ TOV \ We note, however, that the initial 
data used by Shibata et al. and Anderson et al. are dif- 
ferent. Shibata et al. use quasiequilibrium initial data 
for binaries in nearly circular orbits constructed using 
the conformal thin-sandwich (CTS) method. In contrast, 
Anderson et al. set up the initial data by superposing 
boosted metrics of two spherical neutron stars. Ander- 
son et al.'s initial data results in an orbital eccentricity 
of about 0.2 [27|], whereas the eccentricities of the CTS 
initial data are < 0.015 according to a post-Newtonian 
analysis [30| . We point out that a quasi-circular orbit is 
more realistic because gravitational radiation would have 
circularized the orbit long before the binary separation 
reaches a few NS radii. Also, a NS in a quasiequilibrium 
circular orbit will be distorted by tidal effects. 

Anderson et al. report in 12811 that a NSNS with the 
same EOS and masses as in [271 ] but with different initial 
separation leads to a hypermassive neutron star after the 
merger. Later, the star undergoes 'delayed collapse' to a 
black hole due to gravitational radiation. This finding is 
also surprising. The total rest mass M Q I.OSMq ^ is 
smaller than the maximum mass of a uniformly rotating 
star (the supramassive limit), Mq SUP ^ — 1.15Mq TOV ' ) [3l| . 
Hence the remnant cannot be a hypermassive NS, i.e. a 
NS whose mass exceeds the supramassive limit. As a re- 
sult, the star will be unstable to gravitational collapse 
only if a very large amount of angular momentum is re- 
moved. A priori, the expected outcome is that gravi- 
tational radiation removes only some of the angular mo- 
mentum, the remnant acquires some differential rotation, 
and the star settles down to a stable, stationary, rotating 
configuration. 

To better understand the coalescence of NSNSs and 
the role of magnetic fields, we perform a new series of 
simulations using our DLSS code. In this paper, we con- 
sider three models, using the same CTS initial data as 
in [l[ . Specifically, we study the irrotational binary mod- 
els M1414, M1616 and M1418 in [l]. In models M1414 
and M1616, the NSs are of equal mass. In model M1418, 
the ratio of the rest masses of the two NSs are g=0.855. 

We first repeat the calculations of Shibata, Taniguchui 
& Uryu Q for unmagnetized NSNS mergers. Our results 
agree with those reported in Model M1414 results in 
a dynamically stable, differentially rotating hypermassive 
NS. For models M1616 and M1418, the mergers lead to 
prompt collapse to a black hole. The simulations in [l[ 
are terminated soon after black hole formation because of 
grid stretching. We are able to follow the evolution using 
puncture gauge conditions (see, e.g. [12, [HI) until the 
system settles down to a quasi-equilibrium state. This 
allows us to estimate the disk mass around the black 
hole more accurately, and our results are again consis- 
tent with those estimated in We next consider the 
magnetized cases. We add a poloidal field with strength 
B ~ 10 16 G inside each NS of the three NSNS models and 
follow the evolution. While such interior field strength 
may be larger than the value expected for a typical NS, it 



is comparable to the strength inferred for magnetars [34[ 
and is high enough to demonstrate the dynamical ef- 
fects of a magnetic field, if any. For model M1414, the 
merger again results in a differentially rotating, hyper- 
massive neutron star. We see observable difference in 
the magnetized case after the merger as magnetic fields 
are amplified. For model M1616, the system collapses to 
a black hole after the merger as before, leaving negligible 
amount of material to form a disk. This result is unaf- 
fected by the presence of the magnetic field. For model 
M1418, about 1.5% of rest mass is left to form a disk for 
both magnetized and unmagnetized cases. Gravitational 
waveforms for models M1616 and M1418 show only a 
slight difference in amplitude during the entire simula- 
tions. This is because the remnants quickly collapse to 
a black hole after the merger and magnetic fields do not 
have enough time to amplify and alter the dynamics of 
the fluid substantially. This result is not surprising since 
the ratio of magnetic to gravitational potential energy is 
Em/\W\ ~ 10~ 4 initially, and hence the magnetic fields 
are not expected to have an impact on the dynamics be- 
fore they are amplified. 



For models M1414 and M1418, magnetic fields are ex- 
pected to affect the long-term secular evolution of the 
remnants. For the cases where the remnant is a hyper- 
massive neutron star (M1414) , magnetic fields are crucial 
for driving 'delayed collapse' of the star [l||, UU, and the 
resulting remnant could be a central engine for a short- 
hard GRB [3]. The effect of a magnetic field may be 
diminished whenever the merged hypermassive remnants 
develop a bar [25|, The bar leads to dissipation of 
angular momentum by gravitational radiation and may 
result in delayed collapse on a timescale faster than that 
of magnetic field amplification. However, a bar does not 
develop for model M1414. In general, the development 
of a bar depends on the NS EOS. For cases in which 
the remnant consists of a black hole surrounded by a 
disk (M1418), magnetic fields may produce turbulence in 
the disk via MHD instabilities and may generate ultra- 
relativistic jets [H [H, HE [H, [H, \M, El ■ In this paper, 
we are primarily interested in studying the effect of the 
magnetic field during the late inspiral, merging and the 
early post-merger phases, so we do not follow the long- 
term evolution of the remnants. We have previously stud- 
ied the long-term secular evolution of magnetized hyper- 
massive NSs in [l3|, [lj| and the evolution of magnetized 
disks around rotating black holes in [42j |. 



This paper is organized as follows. In Sec [III we briefly 
summarize the basic equations and their specific imple- 
mentation in our GRMHD scheme. In Sec lIHi we present 
the results of our simulations and compare them with 
those in [l|. We summarize our results in Sec. IIVI and 
comment on future directions. 
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II. FORMULATION 
A. Basic equations and numerical methods 

The formulation and numerical scheme for our 
GRMHD simulations are the same as those reported 
in to which the reader may refer for details. Here 

we briefly summarize the method and introduce our no- 
tation. We adopt geometrized units (G = c = 1) except 
where stated explicitly. 

We use the 3+1 formulation of general relativity and 
decompose the metric into the following form: 



ds 2 = 



2 dt 2 + {dx l + (3 l dt) (dx J + (3 J dt) . ( 1 ) 



The fundamental variables for the metric evolution are 
the spatial three- metric 7^ and extrinsic curvature K{j. 
We adopt the BSSN formalism [H, El] to evolve 7^ and 
Kij. In this formalism, the evolution variables are the 
conformal exponent (f> = ln7/12, the conformal 3-metric 



Hi 



-4<p 



jij, three auxiliary functions T l 



-7^, the 



trace of the extrinsic curvature K, and the tracefree part 



of the conformal extrinsic curvature Ai 



-i,b 



(K i: 



jijK/3). Here, 7 = det(7i,). The full spacetime metric 
9fiv is related to the three-metric 7^ by 
n^Tiy, where the future-directed, timelike unit vector 
normal to the time slice can be written in terms of the 
lapse a and shift (3 l as n M = a _1 (l, —/?*). As for the 
gauge conditions, we adopt an advective "1+log" slicing 
condition for the lapse and a second-order "non-shifting- 
shift" HHH] as in our BHNS simulations 

The BSSN equations are evolved with fourth-order ac- 
curate spatial differencing and upwinding on the shift 
advection terms. We apply Sommerfeld outgoing wave 
boundary conditions to all BSSN fields, as in [2Cj. Our 
code is embedded in the Cactus parallelization frame- 
work [46], whereby our second-order iterated Crank- 
Nicholson time-stepping is managed by the MoL, or 
method of lines, thorn. We use the moving puncture 
technique to handle any black hole that may form after 
the merger of the NSNS. The apparent horizon of the 
black hole is computed using the ahf inderdirect Cac- 
tus thorn [47J • Before an apparent horizon appears, we 
find that adding a Hamiltonian constraint term to the 
evolution equation of (f> as in [48[ leads to smaller con- 
straint violation during the evolution. However, when a 
black hole appears, we remove this term as it sometimes 
leads to unstable evolution. 

The fundamental variables in ideal MHD are the rest- 
mass density po, specific internal energy e, pressure P, 
four- velocity u^, and magnetic field measured by a 
normal observer moving with a 4-velocity (note that 
B^n^ — 0). The ideal MHD condition is written as 
u fl F ,lv — 0, where F^ u is the electromagnetic (Faraday) 
tensor. The tensor F^ v and its dual in the ideal MHD 
approximation are given by 



(3) 



where t^ ua i3 is the Levi-Civita tensor. Here we have in- 
troduced an auxiliary magnetic 4- vector 6 M = / -\/47r, 
where B^ is the magnetic field measured by an observer 
comoving with the fluid and is related to B^ by 



(u) 



{5» v + u»u v )B" 
n\u x 



(4) 



The energy-momentum tensor is written as 

rp rynFluid , rp~EM f r\ 

J F — 1 fJ,f ' 1 flL> > \°) 

where Tj^ uld and T™ denote the fluid and electromag- 
netic contributions to the stress-energy tensor. They are 
given by 



T J I J md = Pohu^Uv + Pg^, 



(6) 



and 



(7) 



where h = 1 + e + P/po is the specific enthalpy, and 
b 2 = b^bn. Hence, the total stress-energy tensor becomes 



Ttiv = {poh + b + I P + —j g^ - b^K . (8) 

In our numerical implementation of the GRMHD and 
magnetic induction equations, we evolve the densitized 
density p*, densitized momentum density Si, densitized 
energy density f , and densitized magnetic field B l . They 
are defined as 



P* = -VlPon^, 

Si = -^fT^n^^i, 
t = ^T pLU n^n 1 ' - p* 

B l = ^B\ 



(9) 

(10) 

(11) 
(12) 



Uabp, 



(2) 



During the evolution, we also need the three- velocity v % = 

The MHD and induction equations are written in con- 
servative form for variables p*, S i: f , and B % and evolved 
using an HRSC scheme. Specifically, we use the mono- 
tonized central (MC) scheme [49] for data reconstruction 
and the HLL (Harten, Lax and van-Leer) scheme [5(| 
to compute the flux. The magnetic field B % has to 
satisfy the "no monopole" constraint diB 1 = 0. We 
adopt the flux- in terp olated constrained transport (flux- 
CT) scheme [HI, [52| to impose this constraint. In this 
scheme, the induction equation is differenced in such a 
way that a second order, corner-centered representation 
of the divergence is preserved as a numerical identity. 
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At each timestep, the hydrodynamic "primitive" vari- 
ables (po,P,v l ) must be computed from the "conserva- 
tive" variables (/?*, f , S^). This is done by numerically 
solving the algebraic equations ©-(HI]) together with an 
EOS P = P(po, e)- In this paper, we adopt a T-law EOS, 
P = (r- l> e, with T = 2. 

As in many hydrodynamic simulations in astrophysics, 
we add a tenuous "atmosphere" that covers the compu- 
tational grid outside the star. The atmospheric rest-mass 
density is set to 10 _10 /9 ma x(0), where p max (0) is the initial 
maximum rest- mass density of the stars. As in [l5j |. we 
apply outer boundary conditions on the primitive vari- 
ables p,P,v l , and B l . Outflow boundary conditions are 
imposed on the hydrodynamic variables, and the mag- 
netic field is linearly extrapolated onto the boundaries. 
Finally, the evolution variables p* , Si, and f are recom- 
puted on the boundary. 

Our GRMHD code (DLSS) has been thoroughly tested 
by passing a robust suite of tests. These tests include 
maintaining stable rotating stars in stationary equilib- 
rium, reproducing the exact Oppenheimer-Snyder solu- 
tion for collapse to a black hole, and reproducing analytic 
solutions for MHD shocks, nonlinear MHD wave propa- 
gation, magnetized Bondi accretion, and MHD waves in- 
duced by linear gravitational waves 0. Our DLSS code 
has also been compared with SS's GRMHD code [l0( by 
performing identical simulations of the evolution of mag- 
netized hypermassive NSs fl3l. [Till and of magnetorota- 
tional collapse of stellar cores [ld^We obtain good agree- 
ment between these two independent codes. Our code 
has also been used to study the collapse of very massive, 
magnetized, rotating stars to black hole I21f . evolution 
of merging BHBH |53[ and BHNS binaries [20||, and the 
evolution of relativistic hydrodynamic matter in the pres- 
ence of puncture black holes [54| . Recently, our code has 
been generalized to incorporate (optically thick) radia- 
tion transport and its feedback on hydrodynamic mat- 
ter m. ' 



B. Initial data 

We adopt the same irrotational, quasi-circular NSNS 
initial data as in These initial data set were gener- 
ated by Taniguchi & Gourgoulhon [56j, l57f by numerically 
solving the constraint equations of general relativity in 
the CTS framework. We consider three models studied 
in M1414, M1616 and M1418. 

All models assume an n—1 polytropic EOS for the neu- 
tron stars: P=np\. The compaction, (M^/R)^, is de- 
fined as the ratio of the ADM (Arnowitt-Deser-Misner) 
mass .M* to the areal radius R of a spherical neutron star 
in isolation. For an n—1 polytropic EOS, the compaction 
uniquely specifies the neutron star. We thus label the 
NSNS models by the compaction of each neutron star. 
Model M1418 means the compactions of the two neutron 
stars are 0.14 and 0.18. Hence the two neutron stars do 
not have the same rest masses. For models M1414 and 



M1616, the two neutron stars are of equal rest mass and 
their compactions are 0.14 (for model M1414) and 0.16 
(for model M1616). It is convenient to rescale all quan- 
tities with respect to n. Since k 1 / 2 has dimensions of 
length, we can define the nondimensional variables (58j 
M—n^ 1 / 2 M, R=n~ 1 / 2 R, and p =np . Here M is the 
ADM mass of the binary. The relationship between these 
nondimensional variables and quantities in cgs units are 



M = 10M e 
R = 14.8km 



1/2 



Po 



1.455 x 10 5 cgs 



1.455 x 10 5 cgs 



6.18 x 10 15 gcm 



M 



1/2 



R 



1.455 x 10 5 cgs 



(13) 
(14) 
Po , (15) 



where the value k=1.455 x 10 5 cgs is used by 59]. The 
maximum rest mass for a spherical neutron star for this 

(TOV)_ 



E_OS is ' =0.180, and the maximum ADM mass is 

M™ ax =0.164. Table U summarizes the characteristics of 
our models. 

To study the effects of magnetic fields, we add a small 
poloidal magnetic field to the quasi-equilibrium model. 
We orient our coordinates so that the initial maximum 
densities of the two neutron stars are located at (xi, 0, 0) 
and (x2, 0, 0) with X\ < and Xi > 0. For each neutron 
star I, we specify the magnetic vector potential 



4° = -G/M 2 )4° , 

A« = 0, (J = 1,2) 



A« = (x/w?)A$> , 



(16) 
(17) 



where wi = y/(x — x;) 2 



and we set 



4? 



A b wf 



1 



Po 

.max 

Po 



max(P-P cut ,0) . (18) 



Here A^, n p and P cut are free parameters. The mag- 
netic field in the Z-th neutron star is then computed by 
— n^ l ^ k djA i f > . This guarantees that the magnetic 

constraint djB^ = is automatically satisfied. The pa- 
rameter At, determines the strength of the B-field. The 
cutoff pressure parameter P cut confines the B-field inside 
the neutron star to reside within P > P cu t- The pa- 
rameter n p shifts the location of the maximum B-field in 
the star. A larger value of n p results in the maximum 
B-ficld located in the lower density region of the star. 
In this paper, we set P cut to be 4% or 0.1% of the maxi- 
mum pressure, n p =0 or 3, and set Ab so that the volume- 
averaged magnetic field is 10 16 G(Afo/2.8M Q ), where Mq 
is the total rest mass of the stars. There is no initial ex- 
terior magnetic field. Figure [T] shows the initial magnetic 
field configurations for a widely separated NS compan- 
ion with compaction M*/i2=0.14. The magnetic field 
lines, as well as density distribution, for each NS in our 
binary at t = are slightly distorted from those shown 
in Fig. Q] due to tidal effects. Table [IT] lists the simula- 
tions performed in this paper with a brief summary of 
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TABLE I: Irrotational, quasiequilibrium NSNS models in circular orbit. Here (M*/R)oo is the neutron star compaction, Pq 13 " 
is the maximum nondimensional rest-mass density of a neutron star, Mo is the nondimensional total rest mass of the binary, 
M is the nondimensional ADM mass of the system, J is the ADM angular momentum, q = Mq^/Mq 2 ' is the ratio of the rest 
masses of the stars, Qo is the quasi-circular orbital angular velocity, and Pq — 2tv /fio is the orbital period. 

Model (NU/R)^ pg lax Mo M J/M 2 g MQ Q P /M 

M1414 0.14, 0.14 0.118, 0.118 0.292 0.269 0.951 1.00 0.0326 193 

M1616 0.16, 0.16 0.151, 0.151 0.320 0.292 0.914 1.00 0.0395 158 

M1418 0.14, 0.18 0.118, 0.195 0.317 0.290 0.933 0.855 0.0345 182 



TABLE II: Parameters and results of various runs. The label HNS stands for "hypermassive neutron star," BH stands for 
"black hole," M 1S is the rest mass of the disk around the black hole, and Jh/M^ is the spin parameter of the black hole. 



Run 


Model 


B- field 


n p 


-feu t /-fm ax 


Result 




Jh/M 2 h 


M1414B0 


M1414 


no 






HNS 






M1414B1 




yes 





0.04 


HNS 






M1616B0 


M1616 


no 






BH 


< 10~ 6 


«0.85 


M1616B1 




yes 





0.04 


BH 


< 10~ 4 


^0.85 


M1616B2 




yes 


3 


0.001 


BH 


< 2 x 10~ 4 


^0.85 


M1418B0 


M1418 


no 






BH 


« 0.013 


»0.8 


M1418B1 




yes 





0.04 


BH 


< 0.018 


«0.8 



the outcomes. We note that Anderson et al. [28| also 
set up the magnetic fields in a very similar way. Their 
setup corresponds to setting the parameters n p =0, P cu t 
to the atmosphere value, and Af, such that the maximum 
B-field strength in each NS is 9.6x lO 15 G(1.78Af /M). 

We note that the interior magnetic field strength of 
10 16 G is large compared to values inferred for the surface 
of a typical pulsar (B ~ 10 12 G), but it is comparable to 
the strength in a magnetar [34[. We find that for this 
B-field strength, the magnetic pressure, P mag = b 2 /2, is 
about 0.1% of the gas pressure, and the total magnetic 
energy is ~ 1CP 5 of the ADM mass. We note that the ac- 
curacy of the ADM mass of our CTS initial data is also of 
order 10 -5 [56]. Hence, adding a B-field of this strength 
induces negligible constraint violation and causes only a 
small perturbation to the equilibrium of the stars and 
their orbit. In addition, the magnetic field profile inside 
a neutron star is not known. Our profile is to make the 
ratio P m ag/-P small initially in most regions. As a result, 
the magnetic field introduces only a small perturbation 
to the fluid and no "magnetic wind" is generated in the 
low-density regions of the star. 

In nature, magnetic fields are not confined inside the 
NS, but extend out to the NS exterior. The exterior fields 
of the two NSs in the binary will interact and modify the 
dynamics during the inspiral phase. This effect has been 
studied analyticall y in Newtonian and post-Newtonian 
calculations in [60Ll6ll [62j]. As a rough estimate, we ap- 
proximate the NS by a sphere of radius R. Consider a 
pure dipole exterior field and a nearly uniform interior 
field alligned with the orbital angular momentum. The 
magnetic dipole moment /i is related to the interior field 
strength B by /j, ss BR 3 jl. The accumulated gravita- 
tional wave cycles during the inspiral phase from grav- 



itational wave frquencies / m i„ to / m „ due to magnetic 
dipole interaction is estimated to be [60| 



SN, 



mag 



25 BjBiRfRl 
647 77 2 M 4 



(ttM/) 



-1/3 



(19) 



where rj = M\MijiM\ + M2) 2 is the symmetric mass ra- 
tio of the two NSs. To estimate the upper bound starting 
from our initial data, we set / max = 00 and / m ; n to twice 
of the initial orbital frequency of the binaries in our mod- 
els. We find <5-/V mag < 0.02 for all three models. Hence 
the effect of magnetic dipole interaction during the inspi- 
ral phase is negligible, even assuming our large adopted 
field strength. For this reason, any appreciable dynam- 
ical effects of the magnetic fields will occur only during 
and after the merger phase, for which our confined field 
model will be adequate. 



C. Grid setup 



Even though the dynamics of the system is mainly con- 
centrated in the central region with radius r ^ 10M, we 
set our computational grid to r w 50M in order to ex- 
tract gravitational radiation. To reduce computational 
resources, we employ fisheye coordinates [63, 64] to allo- 
cate the grid more effectively. Fisheye coordinates x l are 
related to the original coordinates x % through the follow- 
ing transformation: 



x 

—r(r), 
r 



r(r) 



ar - 



(1 - a)s cosh[(Y + r )/s] 



■In 



2 tanh(ro /s) cosh[(r — ro)/s] ' 



(20) 
(21) 
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X/R 

FIG. 1: Magnetic field configurations in a widely separated, 
spherical NS companion constructed from the vector poten- 
tial in Eq. (fl8)l for n p =0 and P C ut=0.04P ma x (upper panel), 
and for n p =3 and P cut =0.00f P max (lower panel). Here P max 
is the maximum pressure. The compaction of this NS is 
M»/P=0.14. Dotted (black) concentric circles are rest-mass 
density contours drawn for po/Po mx = 0.9, 0.7, 0.5, 0.3, 0.1 
and 0.001. Solid (green) lines are contours of A v , which co- 
incide with the magnetic field lines in axisymmetry. Contour 
levels of A v are drawn for A v = (A™ ax - Af n )(i/1Q) 2 + A™ 
with i=l, 2, ... 9, where ^™ ax and A™ ln are the maximum 
and minimum value of A v , respectively. 



where r — \/ x 2 + y 2 + z 2 , f = ^Jx 2 + y 2 + z 2 . The 
quantities a, f and s are constant parameters, which 
are set to a = 3, K _1 / 2 f =2.4, and k~ 1//2 s=0.6 in all of 
our simulations. With this choice, the grid spacing in the 
region n~ x l 2 r > 2.4 is increased by a factor of 3. 

We use a cell-centered grid with size 2N x 2N x N in 
x-y-z (assuming equatorial symmetry), covering a com- 
putational domain x e (-NA, NA), y e (-NA,NA), 
and z € (0,iVA). Here N is an integer, A is the grid 
spacing and z is the rotation axis. We set N = 150, 
k _1 / 2 A = 0.04 for models M1414 and M1616, and 
N = 200, k~ 1 / 2 A = 0.03 for models M1418. These val- 
ues of A are chosen so that the diameter of each neu- 
tron star in the equatorial plane is covered by > 40 grid 
points. This is the same resolution used in [lj. We have 
performed a simulation with 75% of the grid spacing (but 
with closer outer boundary) for the unmagnetized M1414 
case and find that the result is very close to our standard 
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t / P 

FIG. 2: Evolution of maximum density pQ lax and minimum 
lapse a m i n for unmagnetized (black solid line) and magne- 
tized (blue dash line) runs of model M1414. The time t 
is normalized by the initial orbital period Po=193M=1.3 x 
10 _5 s (Mo/2.8Mq). The initial maximum rest-mass den- 
sity is pJT x (0)=0.118/k;=7.9 x 10 14 g cm~ 3 (2.8M /M o ) 2 . The 
merger occurs at t ~ lPo. 

resolution. 



III. RESULTS 
A. Model M1414 

Model M1414 is an equal-mass NSNS binary. The total 
rest mass of the system is Mo — 1.6Mq OV . In the ab- 
sence of magnetic field (run M1414B0), the neutron stars 
merge after about one orbit (w 190Af), consistent with 
the result in [l[. We find that magnetic field does not 
change the result. After the merger, the star becomes a 
hypermassive NS. Figure [2] shows the evolution of max- 
imum density p™ ax and minimum lapse a m i n for both 
unmagnetized and magnetized cases. Figure [3] shows the 
density profile along the s-axis and y-axis in the equa- 
torial plane at three different times. We see that the 
unmagnetized case is very close to the simulation in [l| 
(see their Fig. 6a and Fig. 7b). 

Figures 2] and O show the density contours and velocity 
field in the equatorial plane. We see that magnetic field 
causes some mass-shedding in the low-density region. Af- 
ter the merger, we see double cores rotating around the 
center, as in [l[. The star is also pulsating with a large 
amplitude (see Figs. [2] and [3]). These motions give rise 
to gravitational wave signals after the merger (see be- 
low). We note that the apparent bar-like structure seen 
in Figs. S] and [5] in some of the density contours at late 
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FIG. 3: Nondimensional rest-mass density po = f^po along the s-axis (black solid lines) and j/-axis (blue dash lines) in the 
equatorial plane at three different times for the unmagnetized (left) and magnetized (right) cases. 



times is a coordinate (gauge) artifact. In our coordinates, 
the bar is nonrotating and stationary for several periods. 
Figure [5] shows the contours of constant geodesic proper 
distance from the center of the grid. (The geodesic proper 
distance between P and Q in a spatial slice is denned as 
the proper length of the (spatial) geodesic joining P and 
Q). We have verified that the bar- like density contours 
roughly coincide with the constant geodesic proper dis- 
tance contours, indicating that the density distribution 
in the bar-like region is actually close to axisymmetric. 

Magnetic fields do not affect the dynamics of the sys- 
tem prior to merger, as expected. After the merger, 
we see a larger amplitude of pulsation in the mag- 
netized case. However, we expect that the main ef- 
fects of the magnetic field occur on a longer (secu- 
lar) timescale when the field is amplified by differen- 
tial rotation. Magnetic winding will occur on an Alfvcn 
timescale tA^R/vA ~10 rotation periods, where Va = 
\/b 2 /(poh + b 2 ) is the Alfven velocity. Additional am- 
plification will be triggered by the magnetorotational in- 
stability (MM) [TH, [HI . This amplification will lead to 
transport of angular momentum and may trigger a 'de- 
layed' collapse [13L [Til l . We do not follow the evolution 
for that long in this paper. We note that the grid resolu- 
tion in our simulation is insufficient to resolve some of the 
MHD instabilities such as the MRI However, our cur- 
rent resolution is adequate to capture the main effects of 
magnetic fields on our simulation timescale (a few oscil- 
lation periods after the merger). We note, however, that 
we have already performed high-resolution simulations in 
axisymmetry to study the long-term secular evolution of 
magnetized hypermassive NS remnants [l3l . [l5j . We con- 



firm that the magnetic fields cause transport of angular 
momentum, resulting in a delayed collapse of the star to 
a black hole surrounded by a hot, massive disk. 

Figure [7] shows the gravitational waveforms for both 
the magnetized and unmagnetized cases. We compute 
the 1 = 2, m = 2, s = — 2 spin- weighted spherical har- 
monics of the Weyl tensor ip4 at three radii 32 M, 38M, 
and 43M. As shown in Fig. [7£i and 03, the computed 
ip 22 at these three radii are hardly distinguishable when 
properly scaled, indicating that the extracted waveforms 
are being measured in the wave zone and are reliable. 
We see again that the waveforms of the magnetized and 
unmagnetized cases show negligible differences before the 
merger. After the merger, the magnetic field significantly 
affects the motion of the fluid, and the waveforms exhibit 
observable differences in both the amplitude and phase. 
Notice that there are still gravitational wave signals after 
the merger. This is mainly caused by the rotation of the 
double cores and the pulsation of the merged remnant. 

Figures [8] and [9] show the L2 norms of the Hamiltonian 
and momentum constraint violations. The L2 norms are 
normalized as in [2(|. We see that the constraint vio- 
lations are less than 5% during the entire simulations 
(4P = 772M). 



B. Model M1616 

Model M16I6 is also an equal-mass NSNS binary. The 
total rest mass of the system is Mq = 1.78Mq TOV \ In the 
absence of magnetic field, the system merges at t « 150M 
(w 0.9Pq) and the star promptly collapses to a black 




FIG. 4: Snapshots of density and velocity field in the equatorial plane from run M1414B0 (unmagnetized run). Density contours 
are drawn for p /po(0)=0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.01, 0.001, and 0.0001. 



hole. An apparent horizon forms at t = 192M. Fig- 
ure [10] shows the evolution of the maximum density p™ ax 
and minimum lapse a m i n ■ Figure [TTJ shows snapshots of 
the equatorial density contours and velocity field. Our 
result again agrees with |1] . The simulation in [l| is ter- 
minated soon after the formation of an apparent horizon 
because of the grid stretching. Hence [l[ can only give an 
estimate of the upper bound of 0.002Mo for the amount 
of material that can form a disk. We are able to use 
the puncture technique to continue the evolution until 
the system settles down to a stationary state. Figure [F4l 
shows the rest mass of the material outside the apparent 
horizon M out . We see that all the material falls into the 



black hole. The small residual value of M out ~ 10 -6 Mo 
at late times is due to the presence of our (artificial) at- 
mosphere. After t > 250M, the system settles down to a 
vacuum rotating Kerr black-hole spacetime. 

We perform two simulations for the magnetized cases 
with different initial magnetic field profiles (see Table [TTJ). 
Run M1616B1 has the same profile as M1414B1. In run 
M1616B2, more magnetic field is placed in the outer lay- 
ers of the neutron stars, and hence it could counteract 
the gravitational pull of the black hole more effectively. 
We see from Fig. [TO] that runs M1414B1 and M1616B2 
are qualitatively the same as run M1616B0 (unmagne- 
tized run). They both collapse promptly to a black 




FIG. 5: Same as Fig. Q]but for run M1414B1 (magnetized run). 



hole. The apparent horizon appears at about the same 
time (t — 192M) in all three cases. Figures and [T3] 
show snapshots of equatorial density contours and ve- 
locity field. In Fig. [TJl we see that Mo&sk/Mo < 1CT 4 
for both M1616B1 and M1616B2. We see that magnetic 
fields cause a substantial delay in the time at which mate- 
rial in the low-density region falls into the black hole. The 
effect is more pronounced for the case M1616B2 when 
more magnetic field is in the low-density region. The 
simulations of M1616B1 and M1616B2 are terminated at 
t w 500Af when constraint violations start to become 
large (see Fig. HU). However, our current results already 
indicate that even in the presence of magnetic fields the 



amount of material outside the black hole is very small. 

Figure [T5] shows the gravitational waveforms for all 
M1616 models. We see that the waveforms for the three 
runs are very close. This is expected because magnetic 
fields can affect the dynamics substantially only well af- 
ter the merger. However, the merged remnants quickly 
collapse to black holes before the magnetic fields have 
enough time to change the fluid's motion significantly. 
We do see, however, that the amplitude of the waves 
during the collapse are slightly larger for the magnetized 
cases. Figure [T51 shows the two polarizations h + and h x 
of run M1616B0 as observed in the direction 45° to the 
z-axis axis. 

Figures [T71 and IT51 show the L2 norms of the constraint 



10 





FIG. 6: Constant geodesic proper distance from the center, D, 
in the equatorial plane at t — 4Po for run M1414B0. Contours 
are drawn for D/M=1.86i (i=l, 2, 3, 4, 5). 



FIG. 8: Constraint violations for the unmagnetized run of 
model M1414. Upper panel: Normalized L2 norm of the 
Hamiltonian constraint. Lower panel: Normalized L2 norm 
of the x (black solid line) , y (blue dotted line) and z (red dash 
line) components of the momentum constraint. 
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FIG. 7: Gravitation waveforms for model M1414. (a) 
Re(rM^ 22 ) extracted at r = 43M (black solid line), r = 38M 
(blue dotted line), and r = 32M (red dash line), (b) Same as 
(a) but for Im(rMipi 2 ). (c) Re(rMip% 2 ) extracted at r = 43M 
for the unmagnetized (black solid line) and magnetized (cyan 
dash line) cases, (d) Same as (c) but for lm(rMip 22 ). Note 
that in (a) and (b), the lines are hardly distinguishable, show- 
ing good agreement of waveforms at various extraction radii. 




FIG. 9: Same as Fig. |8]but for the magnetized run of model 
M1414. 



violations for runs M1616B0 and M1616B1. The plots 
for M1616B2 are similar to those of M1616B1 and so are 
not shown here. The peaks at t = 192M are due to 
the formation of a central singularity and are contained 
inside the event horizon. After the apparent horizon ap- 
pears at t — 192Af , the constraints are computed only in 
the region outside the apparent horizon, and we see the 



11 




40 80 120 160 
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FIG. 10: Evolution of maximum density pQ lax and minimum 
lapse a min for runs M1616B0 (black solid line), M1616B1 
(blue dotted line) and M1616B2 (red dash line). The merger 
occurs at t ~ 150Af, and an apparent horizon appears at 
t = 192M for both magnetized and unmagnetized cases. 



constraint violations drop to much lower values. This re- 
sult confirms that the large constraint violations in the 
strong-field region are trapped inside the event horizon. 
For run M1616B0, the constraint violations are less than 
3% most of the time. For run M1616B1, the momen- 
tum constraint violations are less than 5% most of the 
time. The Hamiltonian constraint violation is less than 
5% before the collapse, around 5%-10% after the ap- 
parent horizon forms, and gradually increases to 12% at 
t = 500M, after which the simulation is terminated. 



C. Model M1418 

For model M1418, the ratio of the rest masses of the 
neutron stars are q = M (1) /M^ = 0.855. The total 

rest mass of the system is Mq = 1.76Mq° , which is 
about the same as model M1616. As in model M1414, we 
perform an unmagnetized (run M1418B0) and a magne- 
tized (run M1418B1) simulation. The merger occurs at 
t w 180M w IPq. The merged remnant collapses to a 
black hole. An apparent horizon forms at t = 232M for 
both cases. Figure [TH shows the evolution of the maxi- 
mum density p™ ax and minimum lapse a m i n . Figures [2"01 
and !21l show snapshots of equatorial density contours and 
the velocity vector held. 

Figure [22] shows the rest mass of the material outside 
the apparent horizon M out . We see that for the unmag- 
netized case, M out /Mo is settling down to an equilibrium 
value w 0.013, which is consistent with the upper bound 



0.04 reported in [l[. The simulation is terminated at 
t 450M since the constraint violations outside the ap- 
parent horizon increase to more than 15% and so the evo- 
lution becomes inaccurate. We suspect that the problem 
can be solved by increasing the grid resolution near the 
central singularity. We plan to investigate this issue in 
the near future. For the magnetized case, M out /Mo drops 
to 0.018 at the end of our simulation. The simulation is 
terminated at t w 380M when the constraint violations 
exceed 15%. The magnetized run is qualitatively very 
similar to unmagnetized run at this stage. Unlike model 
M1616, however, there is a substantial amount of ma- 
terial left to form a disk in this case. Magnetic fields 
are expected to play an important role in the subsequent 
secular evolution of the disk. The field could drive MHD 
turbulence and generate ultra-relativistic jets. However, 
the disk mass is probably not large enough to produce 
a short-hard GRB in this case [65|]. We have already 
studied axisymmetric, magnetized disk evolution around 
black holes in (42j. 

Figure [23] shows the gravitational waveforms. We see 
that the waveforms of the unmagnetized and magnetized 
simulations are very close, as we have found for model 
M1616. 



IV. SUMMARY 

We have performed a series of new simulations involv- 
ing the coalescence of NSNS binaries in full general rela- 
tivity, using the CTS quasiequilibrium NSNS initial data. 
We considered three models M1414, M1616 and M1418 
previously studied by Shibata, Taniguchi & Uryii [j| . We 
performed unmagnetized and magnetized simulations for 
each model. 

We find that our results for the unmagnetized runs 
agree with those in Q. In particular, the remnant of 
the merger is a hypermassive neutron star for model 
M1414, a rotating black hole with negligible disk for 
model M1616, and a black hole surrounded by a disk 
with a rest mass less than 2% of the total rest mass 
of the system for model M1418. Given our good agree- 
ment with the results of Shibata, Taniguchi & Uryu [l[ 
for un mag netized binaries, the results of Anderson et 
al. [27L |28| ] remain somewhat puzzling for such configu- 
rations. Our simulations on magnetized NSNSs indicate 
that the magnetic fields can cause observable differences 
in the dynamics and gravitational waveforms after the 
merger (especially for model M1414). However, we find 
that the magnetic effects are less significant than those 
reported in [281 ]. 

For model M1414, the merged remnant consists of a 
double core rotating around the center, and the star pul- 
sates. The motion in this remnant emits gravitational ra- 
diation. We see observable difference between the magne- 
tized and unmagnetized cases in the amplitude of the pul- 
sations. Gravitational waveforms also show differences in 
amplitude and phase after the merger. We expect that 
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FIG. 11: Same as Fig. [4] but for run M1616B0 (unmagnetized run). The black region near the center in the last three panels 
denotes the apparent horizon. 



the magnetic field will be amplified by differential rota- 
tion on a (secular) timescale of ~ 10 rotation periods. 

For the more massive model M1616, prompt collapse 
to a black hole occurs following the merger. For the 
unmagnetized case, all the material falls into the black 
hole. For the magnetized cases, we try two different ini- 
tial magnetic-field profiles and find that the final result is 
about the same as the unmagnetized case. We find that 
< 1CP 4 of the total rest mass resides outside the black 
hole at the end of our simulations (500M). We see only a 
slight difference in the amplitude of the gravitational ra- 
diation between the magnetized and unmagnetized cases 
during the collapse. 



For model M1418 consisting of unequal masses, the 
merged remnant also collapses promptly to a black hole 
after the merger, but there is a substantial amount of 
material left to form a disk. We find that the disk mass 
is < 0.013Mo for the unmagnetized case and < O.OI8M0 
for the magnetized case, where Mq is the total rest mass 
of the system. Magnetic fields are crucial for the sub- 
sequent, secular evolution of the disk. We have previ- 
ously performed long-term, axisymmetric simulations of 
magnetized disks around black holes resulting from the 
collapse of hypermassive neutron stars [12] • We find that 
the magnetic fields can cause outflows, depending on the 
EOS and the magnetic field configuration. 
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FIG. 12: Same as Fig, fill but for run M1616B1 (magnetized run). 



In summary, we find that the effects of magnetic fields 
during and shortly after the merger phase are significant 
but not dramatic. We believe that the most important 
role of magnetic fields are on the long-term, secular evo- 
lution of the merged remnants consisting of a hypermas- 
sive neutron star or a black hole surrounded by a disk, 
as we have demonstrated in our axisymmetric simula- 
tions in previous publications [HI, [l5|, [42J. This is not 
to say that a long-term evolution in 3+1 dimensions is 
not necessary. For one thing, a 3+1 simulation of the 
remnants will evolve self-consistently from the NSNS ini- 
tial data. For the other, MHD turbulence is expected 
to be more prominent in three dimensions than in ax- 
isymmetry 66], and hence may affect the dynamics in 



the long-term evolution. We therefore plan to follow the 
long-term evolution of the remnants in 3+1 dimensions 
in the future. We thus need to overcome the difficulty 
of growing constraint violations we observe in our black 
hole evolutions. We suspect the constraint violations can 
be controlled by increasing the spatial resolution near 
the singularity. We currently use a resolution of about 
M/7 near the puncture in our simulations, which is a 
much lower resolution than the resolutions used in most 
recent binary black hole simulations. Higher resolution 
is also required in the case of model M1414 to resolve the 
MHD instabilities. We also plan to perform simulations 
of NSNS binaries with larger initial separation. This will 
allow us to compute the gravitational waveforms with 
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more cycles so that they can be matched with the post- 
Newtonian waveforms. 
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FIG. 22: Rest mass of the material outside the apparent hori- 
zon M out for the unmagnetized (black solid line) and magne- 
tized (blue dash line) runs of model M1418. 
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FIG. 23: Gravitation waveform iff(t — r) for unmagnetized 
(black solid line) and magnetized (blue solid line) runs of 
model M1418. Gravitation waves are extracted at radius 
r = 43M. 



